DMRG and periodic boundary conditions: a quantum information perspective 
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We introduce a picture to analyze the density matrix renormalization group (DMRG) numerical 
method from a quantum information perspective. This leads us to introduce some modifications for 
problems with periodic boundary conditions in which the results are dramatically improved. The 
picture also explains some features of the method in terms of entanglement and teleportation. 
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The discovery and development of the DMRG method 
P, to treat quantum many-body systems has enabled 
us to analyze and understand the physical properties 
of certain condensed matter systems with unprecedent 
precision Originally envisioned for ID systems with 
short-range interactions at zero temperatures, during the 
last years this method has been successfully extended to 
other situations yL Its mathematical foundations have 
been established 0, Q in terms of the so-called matrix 
product states (MPS) and by now there exists a co- 
herent theoretical picture of DMRG. 

At the same time, the field of Quantum Information 
Theory (QIT) has emerged to describe the properties of 
quantum many-body systems from a different point of 
view. A theory of entanglement has been established, 
and has allowed us to describe and understand phenom- 
ena like teleportation 0, and to use them in the fields 
of communication and computation [8| . Recently it has 
been shown that QIT may also shed some new light in 
our understanding of condensed matter s yste ms 0, El , 
and, in particular, in the DMRG method [l2| . 

In this work we analyze the standard DMRG method 
using a physical picture which underlies QIT concepts. 
The picture has its roots in the AKLT model ^3 and 
allows us to understand why DMRG offers much poorer 
results for problems with periodic boundary conditions 
(PBC) than for those with open boundary conditions 
(OBC), something which was realized at the origin of 
DMRG 2] . It also gives a natural way of improving the 
method for problems with PBC, in which several orders 
of magnitude in accuracy can be gained. The impor- 
tance of this result lies in the fact that physically PBC 
are strongly preferable over OBC as boundary effects 
are eliminated and finite size extrapolations can be per- 
formed for much smaller system sizes. 

Let us start by reviewing the simplest version of the 
DMRG method for 1-D spin chains with OBC, which is 
typically represented as B • B 0, 0|- We denote by 
d the dimension of the Hilbert space corresponding to 
each spin, and by D the number of states kept by the 
DMRG method. We assume that the spins at the edges 
have dimension do > D 14]. At some particular step the 
chain is split into two blocks and one spin in between. 
The left block (L) contains spins 1, . . . , M — 1, and the 



right one (R) spins M + 1, . . . , N. Then a set of D x D 
matrices A s are determined such that the state 
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minimizes the energy. The states \a)L,R are orthonormal, 
and have been obtained in previous steps. They can be 
constructed using the recurrence relations 
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(2) 



a' = l s=l 



where the block L' contains the spins 1, . . . , M — 2. The 
new matrices U^- M ^ ,s are determined from A s and fulfill 



d , 



(3) 



For the blocks consisting of the edge spins alone, the \a) 
are taken as the members of an orthonormal set. 

In order to give a pictorial representation of the above 
procedure we introduce at site M two auxiliary D— level 
systems, om and bu- The corresponding Hilbert spaces 
H a j, are spanned by two orthonormal bases \a) a ,b, re- 
spectively. We take L and ajvf (and also R and b]\t) in 
the (unnormalized) maximally entangled state 
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We can always write | = P t 



M 
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)R,b M , where P, 



M 



maps H a ® Hf, — > Hm , with Hm the space corresponding 
to the M-th spin and [cf. (Q)] 



Pm = 



d D 

EE 

s=l a,P=l 
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In fact, we can proceed in the same way at any other site 
k ^ 1, M, N by defining two auxiliary systems and bk 
and a map Qk defined as in but with the matrices U 
instead of the A. For the edge spins 1 and N we define a 
single auxiliary system b\ and a^v, respectively and define 
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FIG. 1: Schematic picture of the DMRG method for the 
B • B (a) and the B • »B (b) configurations. Horizontal lines 
represent maximally entangled states \<j>), the ellipses and cir- 
cles (squares) the operators Q (P) which map the auxiliary 
system into the physical ones. 



accordingly the operators Qi,jv which now map Hb. a ~ * 
i?i,i\r- Thus, the state \& is then obtained by applying 
the operators Qi . . . , , Pm, ■ ■ ■ Qn to the set of maximally 
entangled states </> between the auxiliary systems bk and 
a k+1 (k = 1,...,N- 1) [see Fig. 1(a)]. 

The DMRG procedure can be now represented as fol- 
lows. At location M, one finds an operator Pm acting on 
the subsystems clm and 6jvf by determining the matrices 
A s . From them, one obtains the operator Qm and goes to 
the next step at location M+l. One proceeds in the same 
way, moving to the right, until one reaches the location 
N. At that point, one starts moving to the left until one 
reaches the location 1 at which point it moves again to 
the right. The procedure is continued until a fixed point 
for the energy is reached, something which always occurs 
since the energy is a monotonically decreasing function 
of the step number. This proves that DMRG with the 
B • B is a variational method which always converges. 

The more standard scenario (B • »B) is represented in 
Fig. 1(b). The operator Pm acts on the auxiliary subsys- 
tems om and and maps H a ® Hb — > Hm ® Hm+i- 
In this picture [for both configurations, Figs. 1 (a,b)] it is 
very clear that the two edge spins are treated on a very 
different footing since they are represented by a single 
auxiliary system which is not entangled to any other. 

In the case of a problem with PBC a slight modification 
of the scheme is used 0- The idea is to still separate 
the system into two blocks and two spins as before but 
now with the configuration B • B». This ensures the 
sparseness of the matrices one has to diagonalize and 
thus it increases the speed of the algorithm '2\. One 
can draw the diagram corresponding to this procedure 
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FIG. 2: Proposed configurations for the case of PBC. One 
may also use two spins instead of one 



in a similar way as in Fig. 1. The important point is 
that still there are always two sites (left most and right 
most of both blocks B) which are treated differently since 
they are represented by a single auxiliary spin which is 
not entangled to any other. In our opinion, this is the 
reason of the poor performance of the DMRG method for 
problems with PBC. 

The method we propose is very clear in terms of this 
picture (Fig. 2). One has to substitute at all sites k the 
spin by two auxiliary systems Ofc and bk of dimension 
D, with bk and ak+i (with a^v+i := &i) in a maximally 
entangled states and find the maps Pk : H a ® H\, — > Hk 
which lead to a state 



I*) = Px®P 2 ...P N -l®P N \ 



(6) 



with the minimal energy This minimization can be per- 
formed in a similar way to the one used in the standard 
DMRG method. Before showing how to do this in prac- 
tice, we derive some formulas in terms of these operators. 
We write 



P* = X>}<V>l fc] |, <¥>l fc] l 



(7) 



Thus, the problem is solved once the states ip (or equiva- 
lently, the matrices B) are determined. Note that start- 
ing from these states, it is possible to calculate expecta- 
tion values of products of local observables Q , since 
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where 
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Thus, the main idea to perform the minimization is very 
simple. Given the Hamiltonian H describing the system, 
one chooses one site M and writes the energy as 



E = 



(10) 



where |^ M 1) = © s |^>|, ) is a vector built by concate- 
nating the tps , and Nm and Hm are d x D 2 hermitian 



3 



square matrices which are built using the vectors 
at k 7^ M. For example, Nm — ® s Nq is a block diag- 
onal matrix with identical blocks No which has matrix 
elements (No)(a,a'),(p,p') = ( N o)( a ,0),(c,',p>), with 

n = e[ m+1] ...e[ n] e^...e[ m - 1] . (11) 

Thus, at this step the operator Pm is found by solving 
the generalized eigenvalue problem 

H M \ipM) = \N M \i>M), (12) 

with A minimum, which in turns gives the energy at this 
step. Then one chooses another site and proceeds in the 
same way until the energy converges. At the end we have 
all the Pfc and can evaluate all expectation values. 

The above method is not very efficient numerically. 
First, the matrix Nq may be ill conditioned. Second, one 
stores many matrices (~ TV 2 ) and performs many matrix 
multiplications (~ TV 2 ) at each step. Now we explain 
how one can make the method much more efficient. 

Let assume that we have a set of spins in a ring. The 
idea is to determine operators P k in a clockwise order 
(first Pi, then P2, until Pjv_i), then improve them fol- 
lowing a counterclockwise ordering (from P/v to P2), then 
again clockwise, until the fixed point is reached. At each 
step, a normalization condition similar to is imposed, 
depending on whether we are in a clockwise or counter- 
clockwise cycle, which makes the matrix Nm well be- 
haved. On the other hand, at each step only the opera- 
tors which are strictly needed in later steps are calculated 
in an efficient way and stored. 

The normalization condition is based on the follow- 
ing fact. Given the state 'J, characterized by matrices 
P, if we substitute B^< s -► B^- S X := U^' s and 
B [M+i],s x -i B [M+i], s ^ where x is a nonsingular ma- 
trix, we obtain the same state. Analogously, we can 
substitute BM> s -» FfiM' 8 := M 1 and B^-^ s -» 
£j[M-i],sy-i_ -yy e cnoose x in the clockwise cycles to 
impose (0 and Y in the counterclockwise ones to impose 

W^M'-WyM.') = 1. (13) 

Thus, at the point of determining the operator Pm, 
|*) = Qi®... Qm-i ® Pm ® Qm+i ■ ® Qiv|0) 8Ar (14) 

where and Qfc are defined as in (J7J but with U and 
V instead of P, respectively. Thus, the operators X and 
y are all of them moved over, such that they are now 
included in those corresponding to Pm- It can be eas- 
ily shown that these conditions on the operator U (V) 
are equivalent to imposing that Pi has the maximally 
entangled state \tj>) as right (left) eigenvector with eigen- 
value 1. This is immediately reflected in the fact that the 
matrix Nm is better behaved, which makes the problem 
numerically stable. 



Let us now illustrate how the procedure works with 
simplest nearest neighbor Hamiltonian crl^crif -1-1 ' , namely 
the Ising Model. Let us assume that we are running 
the optimization of the operators clockwise and that we 
want to determine Pm- So far, in previous steps, apart 
from the matrices U and V, we have stored: (a) For each 
k < M, the following four operators: 

r k := 4 ] E? ] ...Et 2] Et 1] , (15a) 
s k := E£}EW...E?- 2] E?- 1 \ (15b) 
t k := E^EfK..Ef^E^- 1 \ (15c) 

h k := ]T P[ 1] P[ 21 . . . fiM+i] . . . E [ k - 2] E[ k -{\5d) 

n=l 

(b) For each k > M other four similar operators which 
contain products from pl fe l to E^ . With them, one can 
build Hm and No by few matrix multiplications and thus 
determine Pm by solving (JT2J. From it, Qm is deter- 
mined. Then, we construct ru-n,s w-i- i^u-i-i and Hm+i 
starting from Tm, , t m and Hm UM- We continue in 
the same vein, finding four matrices at each step, and 
storing them, until we reach N. Then we start moving 
counterclockwise and start constructing the correspond- 
ing four matrices at each step. Notice that in order to 
construct the matrices Hm and Nq we will have to use 
the stored matrices Ijl5(l which were determined when we 
were moving clockwise. Thus, with this procedure we 
have to store of the order of 47V matrices of dimension 
D 2 (apart from the matrices U, V and the last P's) but 
the number of operations per step is independent of N. 
At the end, when we have reached the fixed point, we 
can determine the expectation value of any operator by 
using © and determining the required matrices using 
© . Note that if the problem has translational symme- 
try, then all these evaluations are even simpler. 

We have applied the above method to the spin 1/2 
Heisenberg chain. We have plotted in Fig. 3 the ener- 
gies obtained as a function of D and compared them 
with those obtained by the standard DMRG method with 
OBC and PBC. From the figure it is clear that the ac- 
curacies we obtain are comparable with those obtained 
with DMRG for problems with OBC but much better 
than for PBC. We have determined the errors by com- 
paring with the exact results [lj] . In the insert of Fig. 3 
we have plotted the local bond strength (S' fe 'S' fe+1 l) as a 
function of k. As expected, the result is independent of 
the position k, as opposed to what occurs with OBC. 

Finally we show that the picture introduced here may 
be valuable to understand the properties of states \& in 
terms of the language and tools developed in the field 
of QIT. First, one can easily see that the entropy of the 
block formed by systems (fco, ^0 + 1, . . • , fci) is bounded 
by 21og 2 (P), as this block is connected to the rest only 
via a ko and b kl +i, and thus the rank of the reduced den- 
sity operator for the block is bounded by the product 
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FIG. 3: (Left) Comparison between DMRG (squares) and 
the new method (circles) for PBC, and N = 28. For reference 
the DMRG results for the Heisenberg chain with OBC 
(triangles) are also shown. (Insert) Variation of the local bond 
strength from the average along the chain, calculated with the 
new method and D = 40. 
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FIG. 4: General states can be expressed in the form @. 



of the dimensions of the corresponding Hilbert spaces. 
Secondly, it allows us to show that any state can be writ- 
ten in the form © (MPS 0, 0) if we choose D = d N 
(actually, D = d\~^ is sufficient). We consider at and 
bk as composed of d-level subsystems, a\,...,a^f and 



b k , respectively, and write 



as a tensor product 

* rt n r\ /-j ^ 

fe+1- 



f>d between 6^ and a 



\ , and contains 



of maximally entangled states 
For k = 2, . . . , N, we choose the operators Pk 

where 1%) is a state for all particles but a 
\4>d) for each pair a k -b k (I > k) and |0) for the rest. The 
action of Pk is to teleport the entangled pairs such that 
at the end one has one entangled pairs between the first 
system and all the rest (Fig. 4), while leaving all the other 
auxiliary particles in |0). Finally, the operator P\ is the 
product of two operators. The first acts on particles a\ 
and transforms 10)^ — > Y$>). The second is l Q i ® (771 1, 
where 771) = 0}^ <g> l^d)®^ -1 . This operator first pre- 



pares the desired state *]/ in particles a\ and then uses 
the available entangled pairs to teleport it to the the rest 
of the particles. 

In summary, we have given a pictorial view of the 
DMRG method and have identified the reason of its poor 
performance for problems with PBC. Our picture imme- 
diately leads to a modified version of the DMRG method 
which dramatically improves the results. This is done at 
the expenses of no longer using sparse matrices, some- 
thing which limits its applications. Nevertheless, we be- 
lieve that the method may allow us to treat problems in 
condensed matter systems which so far have been diffi- 
cult to tackle with the standard DMRG method. In any 
case, the present work illustrates how the developments 
made in QIT during the last years may prove useful in 
other branches of Physics. 

We thank M.A. Martin-Delgado for enlightening dis- 
cussions about DMRG and G. Cabrera for sending the 
extended data of [13. Work supported by the DFG (SFB 
631), the european project and network Quprodis and 
Conquest, and the Kompetenznetzwerk der Bayerischen 
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